;;---------------------------------------------------------------------------
;; plot horizonal and zonal mass concentration of major aerosol species 
;; 
;; . zonal mean 
;; . surface 
;; . on pressure levels 
;;
;; kai zhang, max planck institute for meteorology 
;;  
;;---------------------------------------------------------------------------

load "plot_lib.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


begin 

;;---------------------------------------------------------------------------  
;; plot control
;;  
;; kplot  = 0  aerosol mass (no sulfate)  
;;        = 1  sulfur mass  
;;        = 2  soluble number 
;;        = 3  insoluble number
;;---------------------------------------------------------------------------  

  kplot = 1 
  kplot = 0 
  kplot = 2 
  kplot = 4 
  kplot = 5 
  kplot = 6 

  data_type = "zonal" 
  data_type = "surf" 
  plottype  = "pdf" 
  data_class= "burden" 
  data_class= "mass" 
  data_class= "num" 
  data_class= "aod" 

;; experiment name 

  casa = "E0001B" 
  casb = "HAM2R46" 
  casb = "HAM2R46H" 
 
  ddir = "./data/" 
  pdir = "./" 

  ;load "load.ncl" 

  if(data_class.eq."burden") then
     data_type = "surf"
  end if

  plot_ctrl = True
  plot_comp = True
  plot_diff = True
  plot_obs  = True

  if(kplot.eq.5) then 
     plot_obs  = True
  end if 

  if(kplot.eq.6) then 
     plot_obs  = True
  end if 

  plot_msu = False ;;True
  plot_dms = False ;;True
  plot_so2 = False ;;True
  plot_so4 = False ;;True
  plot_mpm = False ;;True
  plot_mdu = False ;;True
  plot_mss = False ;;True
  plot_mbc = False ;;True
  plot_nus = False ;;True
  plot_nks = False ;;True
  plot_nas = False ;;True
  plot_ncs = False ;;True
  plot_nki = False ;;True
  plot_nai = False ;;True
  plot_nci = False ;;True

  if(kplot.eq.0) then
     plot_mpm = True ;;False ;;True
     plot_mdu = True ;;False ;;True
     plot_mss = True ;;False ;;True
     plot_mbc = True ;;False ;;True
     plot_msu = False ;;True
     plot_dms = False ;;True
     plot_so2 = False ;;True
     plot_so4 = False ;;True
     nrow = 4 
     ncol = 3 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_aero"+"_"+casa+"_"+casb 
     lplot = (/plot_msu,plot_mbc,plot_mpm,plot_mdu,plot_mss,plot_dms,plot_so2,plot_so4/) 
     varname_list = (/"msu","mbc","mpm","mso","mdu","mss","dms","so2","so4"/)
     pind  = (/1,2,4,5/) 
  end if

  if(kplot.eq.1) then 
     plot_mpm = False ;;True 
     plot_mdu = False ;;True 
     plot_mss = False ;;True 
     plot_mbc = False ;;True 
     plot_msu = True 
     plot_dms = True 
     plot_so2 = True 
     plot_so4 = True 
     nrow = 4 
     ncol = 3 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_sulfur"+"_"+casa+"_"+casb 
     lplot = (/plot_msu,plot_mbc,plot_mpm,plot_mdu,plot_mss,plot_dms,plot_so2,plot_so4/) 
     varname_list = (/"msu","mbc","mpm","mso","mdu","mss","dms","so2","so4"/)
     pind  = (/6,7,8,0/) 
  end if 

  if(kplot.eq.2) then
     plot_nus = True
     plot_nks = True
     plot_nas = True
     plot_ncs = True
     plot_nki = False ;;True
     plot_nai = False ;;True
     plot_nci = False ;;True
     nrow = 4 
     ncol = 3 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_soluble"+"_"+casa+"_"+casb
     lplot = (/plot_nus,plot_nks,plot_nas,plot_ncs,plot_nki,plot_nai,plot_nci/) 
     varname_list = (/"nus","nks","nas","ncs","nki","nai","nci"/)
     pind  = (/0,1,2,3/) 
  end if

  if(kplot.eq.3) then
     plot_nus = False ;;True
     plot_nks = False ;;True
     plot_nas = False ;;True
     plot_ncs = False ;;True
     plot_nki = True
     plot_nai = True
     plot_nci = True
     nrow = 3
     ncol = 3 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_insoluble"+"_"+casa+"_"+casb
     lplot = (/plot_nus,plot_nks,plot_nas,plot_ncs,plot_nki,plot_nai,plot_nci/) 
     varname_list = (/"nus","nks","nas","ncs","nki","nai","nci"/)
     pind  = (/4,5,6/) 
  end if

  if(kplot.eq.4) then
     plot_aod = True
     plot_abs = True
     plot_ang = True
     nrow = 3
     ncol = 3 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_"+casa+"_"+casb
     lplot = (/plot_aod,plot_abs,plot_ang/)
     varname_list = (/"TAU_2D","ABS_2D","ANG"/)
     pind  = (/0,1,2/)
  end if

  if(kplot.eq.5) then
     plot_aod = True
     nrow = 2
     ncol = 3 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_"+casa+"_"+casb+"_OBS" 
     lplot = (/plot_aod/)
     varname_list = (/"TAU_2D"/)
     pind  = (/0/)
  end if

  if(kplot.eq.6) then
     plot_aod = True
     nrow = 2
     ncol = 2 
     plotname = pdir+"fig_"+data_type+"_"+data_class+"_"+casa+"_"+casb+"_OBS"
     lplot = (/plot_aod/)
     varname_list = (/"TAU_2D","ABS_2D","ANG"/)
     pind  = (/2/)
  end if

  npr = dimsizes(pind) 


  if(data_type.eq."zonal") then
     stra = "         HAM"
     strb = "         HAM2"
     strc = "         HAM2 - HAM" 
  else 
     stra = "     HAM"
     strb = "     HAM2"
     strc = "     HAM2 - HAM" 
  end if 

  case_str = (/stra,strb,strc/) 

  if(kplot.eq.5) then 

     if(data_type.eq."zonal") then
        strd = "         Composite"
        stre = "         MODIS"
        strf = "         MISR"
     else
        strd = "     Composite"
        stre = "     MODIS"
        strf = "     MISR"
     end if

     delete(case_str) 
     case_str = (/stra,strb,strc,strd,stre,strf/) 
  end if 


  if(kplot.eq.6) then
     strd = "     MODIS"
     delete(case_str)
     case_str = (/stra,strb,strc,strd/)
  end if


;; experiment year 

  year = "2000" 


  fna = ddir+"input_"+casa+"_"+casb+"_"+data_class+"_"+data_type+".nc"
  fnb = "nothing" 
  fnc = "nothing" 
  fnd = "nothing" 
  fne = "nothing" 
  fnf = "nothing" 

  if(kplot.eq.5) then  
     fnb = ddir+"OBS_aod_composite_T63_ANN.nc" 
     fnc = ddir+"OBS_aod_sat-2000_T63_ANN.nc" 
     fnd = ddir+"OBS_misr_T63_ANN.nc" 
    ;fnd = ddir+"OBS_aod_aerocom-median_T63_ANN.nc" 
     flb = addfile(fnb,"r")
     flc = addfile(fnc,"r")
     fld = addfile(fnd,"r")
  end if 

  if(kplot.eq.6) then  
     fnb = ddir+"OBS_ang_modis_T63_ANN.nc" 
     flb = addfile(fnb,"r")
  end if


  fla = addfile(fna,"r")
 
  print("% open files: ") 
  print("% ") 
  print("% "+fna) 
  print("% ") 

  vva = fla->ctrl 
  vvb = fla->test
  vvc = fla->diff

  if(kplot.eq.5) then  
     vvd = flb->TAU_2D
     vve = flc->TAU_2D
     vvf = fld->TAU_2D
      vd = vvd(0,:,:)  
      ve = vve(0,:,:)  
      vf = vvf(0,:,:)  
  end if 

  if(kplot.eq.6) then
     vvd = flb->ANG_MODIS
      vd = short2flt(vvd(0,:,:)) 
  end if

  if(data_type.eq."zonal") then 
     va = vva(:,4:,:)  
     vb = vvb(:,4:,:)  
     vc = vvc(:,4:,:)  
  else 
     va = vva 
     vb = vvb 
     vc = vvc 
  end if 

;;----------------------------------------------------------------------
;; convert burden unit from kg/m2 -> mg/m2   
;;----------------------------------------------------------------------

  if(data_class.eq."burden") then 
     va = va * 1.E+6 
     vb = vb * 1.E+6 
     vc = vc * 1.E+6 
  end if 

;;----------------------------------------------------------------------
;; convert number unit from #/kg -> #/cm3 STP 
;;----------------------------------------------------------------------
;; this is done in mkdata_mm.ncl 

;;----------------------------------------------------------------------
;; positive definite 
;;----------------------------------------------------------------------
;; this is done in mkdata_mm.ncl 

  if(data_class.eq."num") then
     va = where(va.le.0.,0.,va)  
     vb = where(vb.le.0.,0.,vb) 
  end if
   
;;----------------------------------------------------------------------
;; open workstation 
;;----------------------------------------------------------------------

  wks   = gsn_open_wks (plottype,plotname)           ; open workstation
 ;gsn_define_colormap(wks,"BlWhRe")               ; choose colormap
 ;gsn_define_colormap(wks,"testcmap")               ; choose colormap
  gsn_define_colormap(wks,"BlueDarkRed18")               ; choose colormap

  np = 18  
  plot = new(np,graphic) 

;; plot settings 

;; font size  

  if(data_type.eq."zonal") then 
     StringFontHeightF  = 0.030
     cnLineLabelFontHeightF = 0.025
     tmXBLabelFontHeightF = 0.030
     tmYLLabelFontHeightF = 0.025
     tiMainFontHeightF  = 0.05
     tiXAxisFontHeightF = 0.035
     tiYAxisFontHeightF = 0.035 ;;35
     lbLabelFontHeightF = 0.025
  end if 

  if(data_type.eq."surf") then
     StringFontHeightF  = 0.035
     cnLineLabelFontHeightF = 0.025
     tmXBLabelFontHeightF = 0.025
     tmYLLabelFontHeightF = 0.025
     tiMainFontHeightF  = 0.05
     tiXAxisFontHeightF = 0.035
     tiYAxisFontHeightF = 0.035 ;;35
     lbLabelFontHeightF = 0.025
  end if
 

  res                      = True                 ; plot mods desired
  res@gsnDraw              = False                ; don't draw
  res@gsnFrame             = False                ; don't advance frame
 ;res@cnLevelSelectionMode = "ManualLevels"       ; manually select levels
 ;res@cnLevelSpacingF      = 5.0                  ; contour spacing
 ;res@cnMinLevelValF       = -50.                 ; min level
 ;res@cnMaxLevelValF       =  50.                 ; max level
  res@cnLineLabelsOn       = False                ; turn on line labels
  res@cnLinesOn            = False                ; turn on line labels
  res@cnFillOn             = True                 ; turn on color fill
  res@gsnSpreadColors      = True                 ; use full range of map
  res@lbLabelAutoStride    = True                 ; optimal labels
  res@tmYRMode             = "Automatic"
  if(data_type.eq."zonal") then 
     res@vpWidthF          = 0.8        ; width of contour plots
  end if 
  if(data_type.eq."surf") then 
     res@mpShapeMode          = "FreeAspect"
     res@vpWidthF             = 1.20
     res@vpHeightF            = 0.60
     res@lbLabelBarOn         = True
  end if 

  res@lbLabelFontHeightF   = lbLabelFontHeightF
 
  if(data_type.eq."zonal") then 
     res@lbOrientation        = "vertical"
     res@pmLabelBarWidthF     = 0.1        ; default is shorter
     res@pmLabelBarHeightF    = 0.6        ; default is taller
  end if 

  if(data_type.eq."surf") then 
     res@lbOrientation        = "horizontal"
     res@pmLabelBarWidthF     = 1.0        ; default is shorter
     res@pmLabelBarHeightF    = 0.1        ; default is taller
     res@pmLabelBarOrthogonalPosF = .15           ; move whole thing down
  end if 

 ;res@cnFillMode           = "RasterFill"  
 ;res@cnRasterSmoothingOn  = True 

;;----------------------------------------------------------------------
;; strings
;;----------------------------------------------------------------------
  res@gsnLeftStringFontHeightF   = StringFontHeightF
  res@gsnCenterStringFontHeightF = StringFontHeightF
  res@gsnRightStringFontHeightF  = StringFontHeightF

;;----------------------------------------------------------------------
;; axis
;;----------------------------------------------------------------------
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF

;;----------------------------------------------------------------------
;; tickmark
;;----------------------------------------------------------------------
  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; start counting 
;;----------------------------------------------------------------------

  ip = 0 

  npv = dimsizes(varname_list) 
  do it = 0, npv-1 
     print("% index : "+it+" varname : "+varname_list(it)) 
  end do 
  print("% ")
  print("% ")

  do is = 0, npr-1 

     id = pind(is) 
     varname = varname_list(id)

     if(plot_ctrl) then
        vp = va(id,:,:)
        plot(ip) = make_plot(wks,vp,res,data_class,data_type,varname,0,case_str,ip)
        ip = ip + 1
     end if

     print("% ctrl ") 

     if(plot_comp) then
        vp = vb(id,:,:)
        plot(ip) = make_plot(wks,vp,res,data_class,data_type,varname,1,case_str,ip)
        ip = ip + 1
     end if

     print("% test ") 

     if(plot_diff) then
        vp = vc(id,:,:)
        plot(ip) = make_plot(wks,vp,res,data_class,data_type,varname,2,case_str,ip)
        ip = ip + 1
     end if

     if(plot_obs) then
        vp = vd(:,:)
        plot(ip) = make_plot(wks,vp,res,data_class,data_type,varname,3,case_str,ip)
        ip = ip + 1
        ;vp = ve(:,:)
        ;plot(ip) = make_plot(wks,vp,res,data_class,data_type,varname,4,case_str,ip)
        ;ip = ip + 1
        ;vp = vf(:,:)
        ;plot(ip) = make_plot(wks,vp,res,data_class,data_type,varname,5,case_str,ip)
        ;ip = ip + 1
     end if

     print("% diff ") 
     print("% ") 
     print("% ") 
     print("% index : "+is+" id : "+id+" "+varname+" "+data_class+" "+data_type)

  end do

;;----------------------------------------------------------------------
;; make panel
;;----------------------------------------------------------------------

  resp                     = True                ; modify the panel plot
 ;resp@txString            = "zonal distribution"
 ;resp@gsnPanelLabelBar    = True                ; add common colorbar
 ;resp@lbLabelFontHeightF  = 0.007               ; make labels smaller
  resp@gsnPanelYWhiteSpacePercent = 5
  resp@gsnPanelXWhiteSpacePercent = 1
  gsn_panel(wks,plot,(/nrow,ncol/),resp)               ; now draw as one plot

  ;system("convert -density 100 "+plotname+"."+plottype+" "+plotname+".png") 
  ;system("convert              "+plotname+"."+plottype+" "+plotname+".pdf") 

  print("% ")
  print("% ")
  print("% ")
  print("% ")

end






